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Abstract 

A  two-dimensional  mechanistic  model  of  a  tubular  solid  oxide  fuel  cell  (SOFC)  considering  momentum,  energy,  mass  and  charge  transport 
is  developed.  The  model  geometry  of  a  single  cell  comprises  an  air-preheating  tube,  air  channel,  fuel  channel,  anode,  cathode  and  electrolyte 
layers.  The  heat  radiation  between  cell  and  air-preheating  tube  is  also  incorporated  into  the  model.  This  allows  the  model  to  predict  heat 
transfer  between  the  cell  and  air-preheating  tube  accurately.  The  model  is  validated  and  shows  good  agreement  with  literature  data.  It  is 
anticipated  that  this  model  can  be  used  to  help  develop  efficient  fuel  cell  designs  and  set  operating  variables  under  practical  conditions.  The 
transport  phenomena  inside  the  cell,  including  gas  flow  behaviour,  temperature,  overpotential,  current  density  and  species  concentration,  are 
analysed  and  discussed  in  detail.  Fuel  and  air  velocities  are  found  to  vary  along  flow  passages  depending  on  the  local  temperature  and  species 
concentrations.  This  model  demonstrates  the  importance  of  incorporating  heat  radiation  into  a  tubular  SOFC  model.  Furthermore,  the  model 
shows  that  the  overall  cell  performance  is  limited  by  O2  diffusion  through  the  thick  porous  cathode  and  points  to  the  development  of  new 
cathode  materials  and  designs  being  important  avenues  to  enhance  cell  performance. 

©  2005  Elsevier  B  .V.  All  rights  reserved. 
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1.  Introduction 

Tubular  solid  oxide  fuel  cells  (SOFCs)  are  considered  to  be 
the  most  advanced  design  for  mid-  to  large-scale  applications 
up  to  2  MW  [1].  However,  due  to  the  complex  interactions  of 
several  physical  phenomena,  rigorous  mathematical  models 
are  required  to  analyse  the  impact  of  cell  design  and  operating 
conditions  on  the  cell  performance. 

Over  the  recent  years,  various  SOFC  models  have  been 
developed.  Campanari  [2]  developed  a  tubular  cell  model 
from  a  thermodynamic  point  of  view  and  included  the  impacts 
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of  thermodynamic  non-idealities  on  cell  performance  in 
terms  of  parameter  correlations.  Nagata  et  al.  [3]  developed 
a  ID  steady- state  tubular  SOFC  model  considering  elec¬ 
trochemical  and  internal  reforming  reactions.  They  divided 
the  cell  into  many  control  volumes  that  behave  as  continu¬ 
ous  stirred  tank  reactors  (CSTRs).  Unfortunately,  the  model 
simplifies  the  actual  cell  geometry  within  a  stack  and  the 
calculation  of  heat  transfer  coefficient  is  also  based  on  fully 
developed  laminar  flow  at  a  constant  heat  flux.  Aguiar  et 
al.  [4,5]  developed  models  for  direct-  and  indirect-internal 
reforming  of  CH4.  However,  they  used  various  correlations 
to  determine  mass  flux,  heat  flux  and  cell  resistance  that 
are  specific  to  particular  operating  conditions.  Moreover, 
they  assumed  uniform  distributions  of  species  concentra¬ 
tions  and  heat  along  the  radial  direction.  This  implies  that  the 
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Nomenclature 

A 

area  perpendicular  to  the  current  flow  (m2) 

B0 

parameter  defined  by  Carman-Kozeny 
equations 

c 

molar  concentration  (kmol  m-3) 

Cp 

molar  specific  heat  (J  kg  mol-1  K-1) 

Dij 

binary  diffusivity  of  a  species  i  into  another 
species  j  (m2  s-1) 

D-v 

1,  k 

effective  Knudsen  diffusivity  of  a  species  i 
through  porous  media  (m2  s-1) 

Di 

mass  diffusion  coefficient  of  species  i  (m2  s-1) 

Deff 

uo2 

effective  diffusivity  through  the  porous 
cathode 

dip 

particle  diameter  (m) 

E° 

reversible  cell  potential  (V) 

F 

Faraday’s  constant  (96  487  A  s  mol-1) 

Fr 

correction  factor  for  heat  radiation 

I An  5  7cat 

current  along  the  peripheral  path  of  anode  and 
cathode,  respectively  (A) 

7e 

current  density  through  each  cell  radial 
segment  (Am-2) 

J 

current  density  (mAcm-2) 

K 

equilibrium  constant 

L 

cell  length  (m) 

^path 

path  length  along  peripheral  segment  of 
electrodes  (m) 

M 

molecular  weight  (kg  kmol-1) 

in 

mass  flux  (kg  m-2  s-1) 

N 

molar  flux  (kmol  m  2  s-1) 

P 

pressure  (Pa) 

Pw 

partial  pressure  of  water  (millibars) 

Q 

volumetric  heating  rate  (W  m-3) 

R 

gas  constant  (8.31441  Jmol-1  K-1) 

r 

radial  coordinate  (m) 

Pc&u  Pa 

n  discretized  electrical  resistance  of  the 
cathode  and  the  anode,  respectively  (Q) 

s 

electrode  interface  conductivity  (£2-1  cm-2) 

Sj 

volumetric  rate  of  generation  of  species  j 
(kmol  m  3  s-1 ) 

T 

temperature  (K) 

U ,  V 

velocity  in  radial  and  axial  directions, 
respectively  (ms-1) 

y 

mole  fraction 

z 

axial  coordinate  (m) 

Greek  letters 

£ 

emissivity  factor 

0An,  0Cat  average  potentials  at  the  anode  and  the 
cathode,  respectively  (V) 

A//eiec 

enthalpy  generation  due  to  electrochemical 
reaction  (kJ  mol-1) 

T] 

overpotential  (V) 

X 

thermal  conductivity  (W  m-1  K-1) 

li  dynamic  viscosity  (N  m-1  s-1) 

0  potential  (V) 

q p  porosity 

p  gas  density  (kg  m- 3 ) 

a  electrical  conductivity  ( Q  - 1  m- 1 ) 

cr0  Stefan-Boltzmann  constant 

(5.67  x  10-8  W m-1  K-4) 

tm,  rK  tortuosity  factors  for  molecular  and  Knudsen 
diffusion,  respectively 

Subscripts  and  superscripts 
AC  air  channel 

An  anode 

APT  air-preheating  tube 

Cat  cathode 

E  electrolyte 

H2  hydrogen 

H2O  water 

N2  nitrogen 

O2  oxygen 


concentration  resistances  caused  by  mass  diffusion  through 
the  porous  electrodes  are  negligible.  The  most  recent  and 
advanced  models  for  a  tubular  SOFC  have  been  presented  by 
Li  and  Chyu  [6]  and  Li  and  Suzuki  [7].  These  models  account 
for  the  actual  geometry  of  the  cell  and  all  major  phenom¬ 
ena,  including  momentum,  heat  and  mass  transport  coupled 
with  electrochemical  and  chemical  reactions.  Unfortunately, 
these  models  neglect  the  diffusion  of  oxygen  through  the 
thick  porous  cathode.  Since  the  O2  diffusivity  is  much  less 
than  the  H2  diffusivity,  the  overall  cell  performance  is  lim¬ 
ited  by  the  O2  diffusion  through  the  porous  electrode  to 
reach  the  reaction  site  at  the  cathode/electrolyte  interface. 
Furthermore,  the  rates  of  the  electrochemical  reactions  are 
expressed  entirely  in  terms  of  exchange  current  densities  as 
adopted  from  Chan  et  al.  [8]  and  are  assumed  to  remain  con¬ 
stant  along  the  cell  length.  This  assumption  can  cause  serious 
errors  in  the  determination  of  the  current  density  profile  along 
the  cell  length  since  the  kinetics  of  electrochemical  reactions 
strongly  depend  on  temperature  and  local  species  concentra¬ 
tions  [3,9,10].  Thus,  the  exchange  current  density  must  vary 
along  the  cell  length.  Finally,  the  effect  of  heat  radiation  was 
not  incorporated  by  Li  and  Suzuki  [7].  Although  the  model 
of  Li  and  Chyu  [6]  considered  heat  radiation,  the  method  for 
incorporating  surface  radiant  exchange  was  unclear.  Campa- 
nari  and  Iora  [11]  developed  a  tubular  SOFC  model  supplied 
by  methane-reformed  gases.  Unfortunately,  the  momentum 
transport  was  not  incorporated  into  this  model  and  the  cal¬ 
culation  of  heat  transfer  coefficients  was  based  on  constant 
Nusselt  number  for  both  the  air  and  fuel  sides,  which  is  valid 
only  when  a  constant  heat  flux  is  transferred  to  the  fluid  flow. 
However,  heat  generation  within  an  SOFC  depends  on  the 
rate  of  electrochemical  reactions  which  is  not  uniform. 
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In  the  present  work,  a  mechanistic  model  of  a  cathode- 
supported  tubular  SOFC  is  developed.  Fuel  and  oxidant  gases 
considered  in  this  model  are  humidified  hydrogen  and  air, 
respectively.  The  model  geometry  comprises  the  air  and 
fuel  channels,  cathode,  electrolyte,  anode  and  air-preheating 
tube.  Momentum,  heat  and  mass  transport  equations  coupled 
with  the  electrochemical  reactions  are  solved  simultaneously. 
Therefore,  correlations  for  heat  and  mass  transfer  are  not 
required.  In  addition,  the  model  incorporates  a  number  of 
potentially  important  factors:  finite  cathode  size,  mass  trans¬ 
port  of  O2  through  the  cathode  using  the  dusty-gas  model, 
dependence  of  the  exchange  current  densities  on  tempera¬ 
ture  and  concentration  and  radiative  heat  transfer  between 
the  cell  and  air-preheating  tube. 


Table  1 

Dimensions  of  cell  components  and  of  the  flow  passages  used  in  the  model 


Item  Value 


I.  Component  thickness  (|xm) 

Supporting  tube  1500 

Cathode  1000 

Electrolyte  50 

Anode  150 

II.  Diameter  (mm) 

Inner  side  of  air-preheating  tube  8.0 

Outer  side  of  air-preheating  tube  9.0 

Inner  side  of  supporting  tube  13.8 

Outer  side  of  anode  19.2 

Outer  diameter  of  fuel  channel  29.2 

III.  Cell  length  (mm)  500 


2.  Model  description 

2.7.  Basic  features  of  the  cell  and  model  geometry 

The  cell  components  used  in  this  simulation  consist  of 
supporting  tube,  cathode,  electrolyte  and  anode  fabricated  in 
a  tubular  shape  with  one  closed  end,  as  illustrated  in  Fig.  1. 
The  cell  dimensions  are  taken  from  the  work  of  Hagiwara  et 
al.  [12]  (Table  1).  For  the  purposes  of  this  model  development, 


energy  and  mass  transport  through  the  supporting  tube  are 
considered  to  be  similar  to  those  through  the  cathode  so  that 
no  distinction  is  made  between  these  components  of  the  cell. 

Due  to  symmetry  conditions,  only  the  portion  on  one  side 
of  the  symmetry  plane  is  considered  in  the  model  (Fig.  2).  The 
anode  thickness  (r-direction)  is  not  taken  into  account  in  the 
model  geometry  since  mass  transfer  and  ohmic  resistances 
through  the  thin  porous  anode  are  much  smaller  than  those 
through  the  other  cell  components. 


Symmetry  plane 


Anode 
Electrolyte 
r 

Cathode 


(a) 


u 


Fuel  inlet 


Air  inlet 


Air-preheating  tube 


Current  out 


Air-preheating  tube 


(b) 


Current  in 


Intereonneet 


Anode 

Electrolyte 

Cathode 

9 


V 


Fuel  channel 


Symmetry  plane 


Intereonneet 


Fig.  1.  Schematic  diagram  of  a  typical  tubular  SOFC:  (a)  longitudinal  section  of  flow  passages  and  (b)  cell  cross-section. 
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r 


4 


Fuel  channel 


Fig.  2.  Schematic  diagram  of  the  model  geometry  (not  to  scale). 


Since  the  variation  of  the  system  properties  in  the 
azimuthal  (0)-direction  is  not  expected  to  be  significant,  the 
model  equations  are  solved  in  two  dimensions  only  in  the  r—z 
plane  (see  Fig.  2)  to  determine  the  profiles  of  current  den¬ 
sities,  temperature  and  species  concentrations  along  the  cell 
length.  However,  since  the  current  flows  in  the  0-direction 
as  it  flows  from  one  interconnect  to  the  other,  an  ohmic  drop 
exists  in  this  direction.  This  effect  is  incorporated  in  the  model 
by  coupling  the  equations  in  the  r-z  plane  to  a  separate  cal¬ 
culation  of  the  ohmic  resistance  along  the  path  of  the  current 
flow  in  the  r-0  plane  at  each  z-position.  Some  simplifications 
and  assumptions  were  made  in  order  to  effectively  link  the 
calculation  in  the  two  planes,  as  described  next. 

2.2.  Model  assumptions 

•  The  model  is  based  on  steady-state,  non-isothermal  oper¬ 
ation  using  humidified  H2  as  the  fuel  source. 

•  Gas  flows  in  the  channels  are  laminar. 

•  The  fuel  cell  operates  with  100%  current  efficiency,  i.e., 
no  reactant  gas  crossover  or  side  reactions. 

•  There  is  no  variation  of  temperature  and  species  concen¬ 
trations  in  the  azimuthal  direction  within  the  cell. 

•  Heat  transfer  within  the  porous  cathode  occurs  by  conduc¬ 
tion  only,  i.e.,  the  enthalpy  carried  by  oxygen  in  the  pore 
space  is  negligible. 

•  Heat  transfer  by  radiation  occurs  only  between  the  air¬ 
preheating  tube  and  cell  structure.  The  radiant  heat  from 
one  cell  to  the  next  cell  within  the  stack  is  negligible 
because  of  geometrical  symmetry,  i.e.,  the  temperature  dis¬ 
tribution  along  the  cell  length  is  the  same  from  one  cell  to 
the  adjacent  cells. 

•  The  mass  transfer  resistance  through  the  porous  anode  is 
ignored  because  of  its  small  thickness  and  comparatively 
high  H2  diffusivity. 

•  The  ohmic  drop  within  the  cell  is  due  to  ionic  resistance 
of  the  electrolyte  and  electronic  resistance  of  the  porous 
electrodes  along  the  peripheral  current  path  in  the  r-0  plane 
at  each  z-position. 

•  The  current  density  along  the  peripheral  path  at  the  cathode 
is  equal  to  that  at  the  anode  since  the  radius  at  the  periphery 
of  the  cathode  is  similar  to  that  of  the  anode. 


2.3.  Equivalent  circuit  of  a  tubular  cell 


In  order  to  describe  the  flow  of  current  and  the  effect  of 
ohmic  resistances  on  cell  performance,  the  equivalent  circuit 
shown  in  Fig.  3  is  used.  The  flow  of  current  from  the  anode 
interconnect  to  the  cathode  interconnect  through  the  cell  com¬ 
ponents  can  be  described  in  terms  of  an  equivalent  circuit  over 
the  cross-sectional  area  of  the  cell.  To  determine  the  ohmic 
drop  within  the  cathode  and  the  anode,  the  peripheral  path  of 
the  cell  is  divided  into  many  small  equal  segments.  The  resis¬ 
tances  of  the  cathode  and  the  anode  to  current  flow  between 
the  nodes  are  given  as: 


L  path,  An 
^An^An 


Fig.  3.  Schematic  diagram  of  equivalent  circuit  to  account  for  ohmic  resis¬ 
tances  in  the  r-6  plane  (not  to  scale). 
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^Cat 


L  path,  Cat 
^Cat^Cat 


where  /path,An  and  /path,cat  are  the  path  lengths  along  each 
peripheral  segment  of  the  anode  and  the  cathode,  respectively. 
o'  An  and  creat  represent  the  conductivities  of  the  anode  and 
the  cathode,  respectively.  However,  in  the  range  of  SOFC 
operating  temperature  (800-1200  °C),  cr An  and  acat  remain 
virtually  constant  [13].  A  An  and  Acat  are  the  areas  perpendic¬ 
ular  to  the  flows  of  I  An  and  /cat,  respectively. 

Assuming  no  ohmic  resistance  within  the  interconnects, 
0a, o  and  0C, a  are  the  same  as  the  potentials  of  the  anode 
interconnects,  i.e.,  0a,o  =  0An  and  0c,A  =  0Cat-  Moreover,  the 
current  density  flowing  from  the  anode  to  the  cathode  (7e)  at 
every  node  is  the  same  because  it  is  determined  by  the  elec¬ 
trode  kinetics  that  are  functions  of  temperature  and  species 
concentrations  for  which  no  angular  dependence  is  assumed. 
The  current  balance  for  the  anode  region  can  be  described 
as: 


For  node  7, y  =  1, . . TV  —  1 : 


0a,  7+1  0a,  j 

R An 


ZeAe  =  0 


For  node  N : 


0a,  A  0a,  A- 1 

^An 


/eAe  =  0 


The  current  balance  of  each  node  at  the  cathode  region  can 
be  formulated  similarly. 


Although  /e  has  no  angular  dependence,  it  varies  in  the 
longitudinal  direction  (z-direction).  It  becomes  computation¬ 
ally  prohibitive  to  solve  the  anode  and  cathode  current  bal¬ 
ances  for  each  value  of  z  and  then  use  them  in  conjunction 
with  the  transport  equations  in  the  r-  and  z-directions.  Con¬ 
sequently,  the  approach  used  is  to  first  solve  the  transport 
equations  in  the  r-  and  z-directions  to  obtain  numerical  values 
for  the  current  density  /An  (=^Cat)  over  the  entire  periph¬ 
eral  path  as  a  function  of  z.  Once  this  is  obtained,  7e  is 
determined  as  the  average  value  of  /ah  over  the  cell  length 
f,  i.e., 

1  fL 

/e  =  —  /  Ja n  d-  (5) 

/  Jo 

The  current  balances  given  by  Eqs.  (3)  and  (4)  are  solved 
once  using  this  value  of  /e  to  yield  one  set  of  anodic  and 
cathodic  potentials  representing  the  average  behaviour  in  the 
axial  direction. 

The  amount  of  current  flowing  through  the  electrolyte 
is  related  to  the  electrochemical  reactions  of  O2  and  H2  at 
the  cathode/electrolyte  and  the  anode/electrolyte  interfaces, 
respectively.  No  homogeneous  reactions  occur  within  the 
electrolyte  and  transport  by  02~  across  the  electrolyte  occurs 
predominantly  by  migration  rather  than  by  diffusion.  Thus, 
the  distribution  of  the  potential,  0e,  within  the  electrolyte  is 


determined  as  follows: 


V  •  (— (TeV0e)  =  0  or 
1  3(ro e  30e /  dr)  9(<7e  #e /  dz) 


This  equation  can  be  solved  numerically  subject  to  the  fol¬ 
lowing  boundary  conditions: 


n  •  (<teV0e)  =  —/An  at  r  =  r e 

(anode/electrolyte  interface)  (7) 


n  •  (crE  V0E)  =  ^Cat  at  r  =  rc 

(cathode/electrolyte  interface)  (8) 


where  /cat  and  /ad  arise  due  to  the  electrochemical  reactions 
of  O2  reduction  and  H2  oxidation,  respectively,  in  units  of 
A  m-2.  These  rate  expressions  are  obtained  from  the  work  of 
Nagata  et  al.  [3]: 


10  0007?/ 
3 ~F 


^An 


—  exp 


100007?/ 
4 F 


^Cat 


f  ~2ACat\ 
\  RT  ) 


—  exp 


(^)l 

(10) 


where  r/An  and  77 cat  represent  the  anode  and  cathode  overpo¬ 
tentials,  respectively,  i.e., 

OAn  —  0An  0E  ^An  (H) 

OCnt  =  0Cat  -  0E  -  ^Cat 

where  0An  and  0cat  are  the  average  potentials  at  the  anode 
and  the  cathode,  respectively,  estimated  by  Simpson’s  rule 
using  the  potentials  obtained  from  the  equivalent  circuit: 


0An  +  40aq  +  20a, 2  +  40a,3  +  •  •  • 


+  20a,7V- 2  +  40a,  A— 1  +  0a, 

3N 


(13) 


0C, 0  +  40c,  1  +  20c, 2  +  40c, 3  +  •  •  • 


+  20c,  N-2  +  40c,  N- 1  +  0Cat 
3N 


(14) 


Since  the  electrode  resistances  remain  essentially  constant 
over  the  range  of  temperatures  within  the  SOFC,  the  poten¬ 
tials  within  the  electrodes  can  be  determined  once  /e  is  spec¬ 
ified.  The  dependence  of  0An  —  0An  and  0cat  —  0Cat  on  7e  at 
temperatures  between  800  and  1200  °C  is  shown  in  Fig.  4. 
These  simple  relationships  were  used  to  facilitate  the  numer¬ 
ical  algorithm  to  solve  the  model. 
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Fig.  4.  Plot  of  0An  —  0An  and  0cat  —  0Cat  as  a  function  of  /g  between  800 
and  1200  °C. 


^An  and  E°Cat  are  the  reversible  half-cell  potentials  of  the 
anodic  and  cathodic  reactions,  respectively,  as  expressed  by: 


/7  0  _ 

^  An  — 


17O  _ 

^Cat  — 


RT 
4 J 

4F 


Ilf  P02,An 


In  J7Q2,Cat 


(15) 


(16) 


where  po2:An  and  j7o2,Cat  are  partial  pressures  of  O2  (atm) 
at  the  anode/electrolyte  and  cathode/electrolyte  interfaces, 
respectively.  The  former  can  be  calculated  assuming  that 
Po2,An  is  in  chemical  equilibrium  with  other  fuel  gases,  and 
hence  po2, An  can  be  written  as: 


Momentum  in  ^-direction: 

d(PiVjVi)  1  djrpiUjVi) 

3  z  r  dr 


—  _i_  JL 

3  z  3  z 


Pi 


3 \h 
3  z 


13/  3  Vf 

+  r  dr  [  nl  ‘  97 


3  /  3uz’  \  13/  3 m 

+  Jz{IXi5/)+~rd~r\riXi57 

Momentum  in  r-direction: 


3  (PiVjUj)  +  1  3  (rpjUjUi) 

dz  r  dr 

_  _  3 p;  3  /  3 M,- 

3r  3z  \  z  3z 


13/  dui 
+  ~rdr  (mi57 


3  /  3n/  \  13/  dui 

+  Tz  Vi5l)  +  -rd~r{mi57 

Energy  transport: 

dipiCpAJ)  1  djrpjCp.UjT) 

dz  r  dr 


fJ  i  U  i 


3/3  T 


dz  V  ' 

Species  transport: 


13/  3T 

T - (  rA/  — 

r  dr  \  dr 


+  G 


3  (cjMjv{)  1  3  (rcjMjUi) 

dz  r  dr 


3 

3z 


D 


3c 


j 


dz 


1  3 

r  dr 


rD 


J 


dcj_ 

dr 


+  S 


j 


(21) 


(22) 


(23) 


(24) 


/?02,An 


J7H2Q,An 

^H2jPH2,An 


(17) 


where /Th2  is  the  equilibrium  constant  for  H2  +  I/2O2  — >  H2O 
which  depends  on  temperature. 

The  terms  scat  and  SAn  given  in  Eqs.  (9)  and  (10)  are 
estimated  by  the  following  correlations  [3]: 


sCat  =  1.047  x  106  exp 

^7 

sAn  =  6.166  x  10  exp 


130  000  N 

\  0.50 

j  P 02,Cat 

(18) 

RT  j 

120  000  Z 

I  0.133 

1  P  02,An 

(19) 

RT  ) 

2.4.  Momentum,  mass  and  energy  transport  equations 

According  to  the  geometry  shown  in  Fig.  2,  the  equations 
for  momentum,  heat  and  mass  transport  as  indicated  in  Eqs. 
(20)-(24)  are  solved  simultaneously  within  the  gas  flow  pas¬ 
sages  (air-preheating  tube,  air  channel  and  fuel  channel). 


where  i  denotes  fuel  or  air  and  j  represents  H2  in  the  fuel 
stream  or  O2  in  the  oxidant  stream.  The  mass  diffusion 
coefficients,  Du2  and  Do2 ,  as  well  as  the  air  and  fuel  heat 
capacities,  are  determined  using  the  equations  found  in  the 
literature  [14,15]. 

The  sink  term  Q  in  Eq.  (23)  and  the  source  term  Sj  in  Eq. 
(24)  are  set  to  zero  because  generation/consumption  of  heat 
and  mass  is  considered  to  occur  only  at  the  domain  bound¬ 
aries.  In  this  model,  O2  and  H2  consumption  takes  place  at 
the  boundaries,  r=rc  and  rg,  respectively.  Heat  generation 
due  to  the  electrochemical  reaction  is  assumed  to  take  place 
at  r-  rg. 

Within  the  porous  cathode,  the  momentum  and  mass  trans¬ 
port  equations  are  different  from  the  previous  equations  since 
the  porous  cathode  is  treated  as  a  barrier  to  gas  transport.  On 
the  other  hand,  the  energy  equation  can  be  applied  inside  the 
porous  cathode  without  much  modification.  Only  the  correc¬ 
tion  of  the  thermal  conductivity  to  combine  the  properties  of 
the  cathode  material,  lanthanum  strontium  manganite  (LSM) 
and  the  gas  filling  the  pores  must  be  made: 


Overall  mass  balance: 


/-Cat  =  ^CatZair  +  (1  “  </Cat)/LSM 


(25) 


d(pjVj) 

dz 


1  3  (rpM) 
r  dr 


(20) 


where  <^cat  is  the  cathode  porosity  and  Zlsm  is  the  thermal 
conductivity  of  the  cathode  material. 
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The  most  practical  model  widely  used  to  describe  mass 
transport  inside  porous  media  is  the  dusty -gas  model  [16-18]. 
This  model  allows  for  both  Knudsen  and  bulk  diffusion  mech¬ 
anisms  in  the  overall  diffusive  flux  and  mass-average  velocity 
in  convective  flux.  The  expression  governing  the  molar  flux 
of  gas  is  given  by: 


~P 

RT 


Nj  ^  ymNi  —  yiNm 

TSff  +  Z_^  Adf 
IX  l=\,mfl  l-m 


+ 


yi 

RT 


(26) 


where  yi  and  p  denote  mole  fraction  of  species  /  and  the  total 
pressure  inside  porous  media,  respectively.  Df[  and  DfJm 
are  the  effective  Knudsen  and  molecular  diffusivities,  respec¬ 
tively.  Ni  and  Nm  are  the  molar  fluxes  of  the  participating 
gases.  B0  characterizes  the  pore  structure  and  can  be  deter¬ 
mined  from  the  Carman-Kozeny  relationship  [19] 


(27) 


where  <ip,  r  and  cp  are  the  particle  size,  tortuosity  factor 
and  the  porosity  of  the  porous  media,  respectively.  The  first 
and  second  terms  on  the  right-hand  side  of  Eq.  (26)  repre¬ 
sent  Knudsen  and  molecular  diffusions,  respectively,  while 
the  third  term  accounts  for  convective  flow  as  described  by 
Darcy’s  law. 

With  typical  pore  sizes  and  material  properties  of  the 
porous  cathode  ( dv  =  1  jxm,  r  =  4,  cp  =  0.3),  B0  is  estimated 
to  be  1.91  x  10-16  and  1  +  B0p/Df^/i  to  be  close  to  1.  This 
indicates  that  convective  transport  inside  the  porous  cathode 
is  insignificant.  Substituting  this  approximation  into  Eq.  (26) 
an  expression  used  to  determine  O2  mass  flux  is  as  follows: 


cathode  porosity.  Iwata  et  al.  [10]  used  the  same  value  of 
tm  and  tk  in  the  range  of  3-6  in  their  mathematical  model. 
Williford  et  al.  [20]  used  a  Wicke-Kallenbach  experimental 
apparatus  to  estimate  the  tortuosity  of  porous  anode.  They 
found  that  the  anode  materials  with  typical  porosities  greater 
than  30%  exhibit  tortuosities  in  the  range  of  2. 5-3.0,  which  is 
somewhat  smaller  than  what  Iwata  et  al.  expected.  However, 
in  the  present  model,  the  tortuosity  factors  are  estimated  from 
formula  based  on  Monte-Carlo  simulations  [21]. 

Since  the  O2  reduction  reaction  occurs  only  at  the  cath¬ 
ode/electrolyte  interface,  the  steady-state  mass  transport  for 
O2  within  the  cathode  layer  yields 

V  •  rho2  =  0  (32) 

2.5.  Boundary  conditions 

On  the  symmetrical  axis  at  the  air-preheating  tube,  i.e., 
r  =  0:  wair  =  0,  3T/dr  =  0,  dco2/dr  =  0  and  dv^/dr  =  0,  and  at 
the  fuel  channel,  i.e.,  r  =  rpc5  Wfuel  =  0,  dT/dr  =  0 ,  dcn2/dr  =  0 
and  dvfUQ\/dr  =  0.  At  the  fuel  inlet,  i.e.,  z  =  0,  the  inlet 
velocity,  temperature  and  species  mole  fraction  are  pre¬ 
scribed  by  the  known  operating  conditions.  Similar  inlet 
conditions  for  air  are  given  at  z  =  L.  At  the  outlet  of  the 
air  and  fuel  channels  (i:  fuel  or  air)  at  z  =  L,  pi  -  1  atm, 
dui/dz  =  0  and  dvi/dz  =  0.  For  heat  and  species  transport, 
3779z  =  0,  dco2/dz  =  0  and  dcu2/dz  =  0.  At  the  closed 
end  (z  —  0),  uSL{v  =  0,  nair  =  0,  T=  7fueiiniet  and  Vco2(6  z)  — 
0.  Since  the  electrochemical  reactions  consume  O2  and 
H2,  the  boundary  conditions  for  O2  and  H2  transport  are 
given  by  n  •  mo2  =  — /cat47o2/4E  at  r-rc  and  n  •  m\\2  = 
—  J\nM\\2/2F  at  r=re,  respectively.  For  transport  within 
the  electrolyte,  the  boundary  conditions  are  n-(crE  V0e)  =  /cat 
and  n-(crEV(/>E)  =  — ^An  at  r=rc  and  r£,  respectively.  It  is 
assumed  that  the  heat  generation  in  the  tubular  cell  is  due 


m  o2  = 


£>eff  £>eff 
^02,k^02-N2 


^02-N2  +  (1  c02/Ctot)7)o2,k  (c02 / ctot ) ^N2 , k 


eff 


»eff 


where  ct0 1  is  the  total  gas  concentration  in  the  porous  cathode 

(ctot  =  cn2  +  co2).  £>oTn2’  Dol,k  and  DN2,k  can  be  deter_ 
mined  as  follows: 


D 


eff 

o2— N2 


=  A> 


DoU  =  Oo2,k  ( 


^N2,k  —  £*N2,k 


(29) 


(30) 


(31) 


Fr  =  t 


M02Vc0: 


(28) 


to  the  exothermic  electrochemical  reaction  at  r  =  r^.  Conse¬ 
quently,  n-(AEV7)  =  /An[A//eiec/2E+(0cat -0An)]  applies 
at  r=r£,  where  A//eiec  is  the  enthalpy  change  due  to  the 
overall  electrochemical  reaction,  i.e., 


A Hq\qC  =  AHyu2o  -  AHf  n2  -  0.5  AHf  o2 


(33) 


The  radiative  heat  exchange  between  the  surface  of  the  cath¬ 
ode  (r=r ac)  and  that  of  the  air-preheating  tube  (r=  tapt) 
is  incorporated  as  boundary  conditions,  n-(A.APTVr)  = 
cr0FY(T4\r=rAr  -  T4|r=rAPT)  at  r=rAc  and  n  •  (AaptVT)  = 


—<70Fy(T4  \r=rAC  —  T4  |r=r  )  at  r  =  tapt 9  where  Fr  denotes 


4 


zrAC 


=^APT 


the  correction  factor  for  heat  radiation  defined  as  [22] : 

^air^APT(l  £air) 


1  1  1  T  (^Cat^air  £air  ^Cat)  f  (1  £air)(l  ^APT) 


(34) 


where  tm  and  tk  denote  the  tortuosity  factors  for  molecular 
and  Knudsen  diffusion,  respectively,  and  are  functions  of  the 


where  £air,  £apt  and  £cat  represent  the  emissivities  of  air, 
air-preheating  tube  and  cathode,  respectively.  £air  strongly 
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depends  on  the  amount  of  H2O  in  air  as  given  by  [23]: 
eair  =  0.66  +  0.4  (35) 

where  Pw  is  the  partial  pressure  of  H2O  (millibars). 

2.6.  Numerical  implementation 

The  finite-element  software  package  FEMLAB®  [24]  was 
used  to  solve  the  transport  equations  coupled  with  electro¬ 
chemical  reactions.  Depending  on  the  initial  guesses  for  the 
state  variables,  the  time  required  for  each  simulation  was  from 
10  to  40  min  on  a  Pentium  IV  2.6  GHz.  It  should  be  noted  here 
that  the  number  of  grid  points  used  to  discretize  the  geometry 
is  a  compromise  between  accuracy  of  model  results  and  com¬ 
putational  capacity.  The  grid  spacing  at  the  anode/electrolyte 
and  cathode/electrolyte  interfaces  along  the  cell  length  was 
finer  than  in  other  regions  since  high  gradients  of  temperature 
and  species  mole  fractions  occur  in  this  region. 

3.  Results  and  discussion 

3.1.  Model  validation 

To  ensure  that  the  model  can  predict  the  observed 
behaviour  of  an  operating  tubular  cell,  a  simulation  was  done 
to  compare  with  the  experimental  data  of  Hagiwara  et  al. 
[12].  Physical  parameters  of  cell  components  and  cell  oper¬ 
ating  conditions  used  in  this  simulation  are  summarized  in 
Tables  2  and  3,  respectively. 

In  this  simulation,  the  values  of  fuel  and  air  flow  rates 
were  adjusted  to  match  the  values  of  utilization  factors  given 
in  Table  3.  For  different  operating  cell  potentials,  they  were 
found  to  vary  from  ~0.002  to  0.008  kgh-1  for  the  fuel  and 
^0.4  to  0.9  kgh-1  for  air.  However,  it  is  impossible  to  com¬ 
pare  these  values  with  the  experimental  data  since  they  were 
not  reported  by  Hagiwara  et  al.  [12].  The  current-potential 
plot  obtained  from  this  simulation  is  compared  with  the 
experimental  data  and  the  predictions  from  the  model  of  Li 
and  Suzuki  [7]  in  Fig.  5.  The  present  model  can  predict  the 
experimentally  observed  dependence  of  current  density  on 


Table  3 


Operating  conditions  used  in  this  simulation 


Item 

Symbol 

Value 

Species  and  its  molar  fraction 

Hydrogen 

cH2,in 

98.6% 

Water 

CH2  O,  in 

1.4% 

Nitrogen 

^N2,m 

79.0% 

Oxygen 

c02,in 

21.0% 

Inlet  pressure 

Fuel 

Pfuel,in 

1.0135  x  105  Pa 

Air 

Pair, in 

1.0135  x  105Pa 

Inlet  temperature 

Fuel 

Tfuel,in 

870.0  °C 

Air 

T, air, in 

600.0  °C 

Fuel  utilization 

P  fuel 

85.0% 

Air  utilization 

C  ai  i‘ 

16.7% 

Cell  potential 

Lcll 

0.6-0. 8  V 

Average  current  density 

1e 

100-520  mA  cm-2 
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Fig.  5.  Comparison  of  current-potential  plots  by  the  present  model  (♦), 
model  from  Li  and  Suzuki  [7]  (•)  and  the  experimental  data  (Hagiwara  et 
al.  [12])  (0)  at  t/air  ~  0.167  and  t/fuei  ~  0.85. 

cell  potential  very  well.  Excellent  agreement  is  also  obtained 
in  the  low  current  density  region  where  the  fit  of  the  model 
from  Li  and  Suzuki  [7]  shows  some  deviation.  As  expected, 
higher  current  density  leads  to  a  lower  cell  potential  due  to 


Table  2 


Physical  properties  of  the  cell  components  used  in  this  simulation 


Item 

Symbol 

Value 

Reference 

Cathode  porosity 

<?Cat 

0.30 

[8] 

Tortuosity  for  molecular  diffusion 

1.50 

[21] 

Tortuosity  for  Knudsen  diffusion 

3.44 

[21] 

Ionic  conductivity  of  electrolyte 

or 

3.3  x  104 exp(— 10300/D m-1 

[13] 

Electronic  conductivity  of  cathode 

^Cat 

12  800  m-1 

[13] 

Electronic  conductivity  of  anode 

&  An 

30  300S2-1  m-1 

[13] 

Thermal  conductivity  of  cathode 

kcat 

4Wm“1r1 

[26] 

Thermal  conductivity  of  electrolyte 

7-e 

2  Wm_1  K-1 

[26] 

Thermal  conductivity  of  anode 

■^An 

2  Wm_1  K-1 

[26] 

Thermal  conductivity  of  preheating  tube 

^APT 

l.lWm-1  K”1 

[6] 

Thermal  emissivity  of  preheating  tube 

£apt 

0.05 

[25] 

Thermal  emissivity  of  cathode 

£Cat 

0.35 

[25] 
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Fig.  6.  Surface  and  arrow  plots  of  velocity  profiles  for  air  (a)  and  fuel  (b)  inside  flow  passages  for  the  base  case  model  (/g  ~  300  mA  cm  2,  Vceii  ~0.7  V, 
C/air  ~  0.167  and  t/fue i  ~  0.85). 


the  increase  in  the  activation,  ohmic  and  concentration  over¬ 
potentials. 

3.2.  Simulation  results 
3.2.1.  Flow  behaviour 

The  velocity  profiles  for  air  and  fuel  are  presented  in  Fig.  6. 
The  gas  velocity  inside  the  air-preheating  tube  increases 
along  the  tube  length  since  air  expands  due  to  preheating 
(Fig.  6(a)).  Once  air  is  discharged  at  the  closed  end  of  the 
cell,  it  impinges  on  the  cell  wall  and  flows  further  along  the 
air  channel  in  the  reverse  direction  in  the  annular  region 
between  the  cell  and  the  air-preheating  tube.  At  a  given 
mass  flow  rate,  the  velocity  inside  a  tube  is  inversely  propor¬ 


tional  to  the  tube  diameter.  Accordingly,  the  lower  velocity 
at  the  closed  end  of  the  cell  is  expected  since  the  size  of 
the  closed  end  tube  is  larger  than  that  of  the  air-preheating 
tube.  Moreover,  once  air  hits  the  cell  wall,  it  loses  momen¬ 
tum  and  its  velocity  is  reduced.  The  velocity  profile  changes 
only  slightly  along  the  cell  length  outside  the  air-preheating 
tube  because  the  temperature  does  not  change  much  in  this 
region.  This  is  confirmed  in  the  temperature  profile  along  the 
cell  length  as  shown  in  Fig.  7.  Although  the  mass  of  O2  is  con¬ 
sumed  by  the  electrochemical  reaction  along  the  cell  length, 
this  does  not  affect  the  air  velocity  profile  since  the  total 
mass  of  air  flowing  through  the  air  channel  is  much  greater 
than  the  mass  of  oxygen  consumed  by  the  electrochemical 
reactions. 


600  700 


800 


900  1000  1100 


Z  (ni) 


Fig.  7.  Surface  and  contour  plots  of  cell  temperature  distribution  (°C)  for  the  base  case  model  (7g  ~  300  mA  cm  2,  Vceii  =  0.7  V,  t/air  ~  0. 167  and  t/fuei  ~  0.85). 
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Fig.  8.  Rate  of  radiative  heat  transfer  from  the  surface  of  the  cell  to  the 
surface  of  the  air-preheating  tube  along  the  cell  length  at  r=  tac- 

In  contrast,  the  velocity  profile  of  the  fuel,  as  shown  in 
Fig.  6(b),  decreases  along  the  cell  length  although  the  tem¬ 
perature  along  the  cell  length  increases.  This  can  be  explained 
by  the  fact  that  H2  is  being  consumed  electrochemically  to 
form  water  in  the  fuel  channel.  Although  the  molar  consump¬ 
tion  of  H2  is  the  same  as  the  molar  generation  of  H2O,  the 
rate  of  H2O  generation  on  a  mass  basis  is  nine  times  the  rate 
of  H2  consumption.  As  a  result,  the  velocity  of  fuel  stream  is 
reduced  because  of  an  increase  in  its  total  mass. 

3.2.2.  Temperature  distribution 

The  temperature  distribution  is  crucial  to  SOFC  opera¬ 
tion.  The  reaction  kinetics  strongly  depends  on  the  local  cell 
temperature,  thus  affecting  the  overall  cell  performance.  Fur¬ 
thermore,  heat  management  is  a  challenging  problem  during 
cell  operation.  Since  the  cell  length  is  quite  long,  cracking 
of  the  ceramic  material  can  occur  due  to  thermal  gradients 
and  hot  spots  within  the  cell.  Consequently,  this  simulation 
focuses  not  only  on  the  temperature  distribution  within  the 


Fig.  9.  Comparison  of  computed  cell  temperature  profiles  along  the  cell 
length  between  the  cases  where  heat  radiation  is  included  and  excluded. 


cell,  but  also  on  improving  cell  design  so  as  to  minimize 
temperature  non-uniformity. 

Fig.  7  shows  the  temperature  distribution  inside  the  cell. 
In  the  air-preheating  tube,  the  air  temperature  increases  as 
one  approaches  the  closed  end  of  the  cell.  The  gradient  is 
very  high  near  the  air  inlet  and  low  at  the  closed  end  of  the 
cell  before  it  enters  the  air  channel.  This  can  be  explained  by 
the  fact  that  a  high  amount  of  heat  is  transferred  from  the  cell 
to  the  preheating  tube  due  to  the  exothermic  electrochemical 
reaction.  A  significant  rate  of  heat  transfer  occurs  near  the 
inlet  of  the  preheating  tube  because  of  a  large  temperature 
difference  between  the  air  and  the  cell.  When  this  difference 
is  reduced,  the  rate  of  heat  transfer  also  declines  as  observed 
near  the  closed  end.  The  temperature  variation  within  the  air 
channel  is  less  than  that  inside  the  air-preheating  tube  since 
the  temperature  difference  between  the  air  channel  and  the 
cell  is  not  as  high  as  that  between  the  air-preheating  tube 
and  the  cell.  The  radial  temperature  distribution  within  the 


z  (m) 

Fig.  10.  Distributions  of  £°  (a),  overpotential  (b)  and  current  density 
(c)  along  the  cell  length  for  the  base  case  model  (/g  ~  300  mA cm-2, 
Kell  =  0.7  V,  Uak  ~  0. 167  and  t/fue  1  ~  0.85). 
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fuel  channel  and  cell  structure  is  quite  uniform  primarily  due 
to  high  thermal  conductivity  of  the  electrode  material.  With 
regard  to  the  axial  direction,  the  temperature  within  the  fuel 
channel  reaches  a  maximum  of  about  ~  1 1 10  °C  near  the  cen¬ 
tre  of  the  cell,  which  is  ^230  °C  greater  than  the  inlet  fuel 
temperature.  Therefore,  this  region  requires  particular  atten¬ 
tion  since  the  electrode  material  can  crack  due  to  the  stress 
created  by  hot  spots. 

The  three  mechanisms  that  drive  heat  from  the  cell  to  the 
air  are  conduction,  convection  and  radiation.  Radiative  effects 
are  significant  only  when  the  operating  temperature  is  quite 
high,  as  in  a  SOFC.  Fig.  8  shows  the  amount  of  heat  radi¬ 
ated  from  the  cell  surface  to  the  surface  of  the  air-preheating 
tube.  Near  the  fuel  inlet,  the  temperature  difference  between 
the  air-preheating  tube  and  the  cell  is  low  enough  that  only  a 
small  amount  of  heat  is  transferred  to  the  air-preheating  tube. 
On  the  other  hand,  it  is  relatively  uniform  between  z  ~  0.26 
and  ~0.34m  where  the  temperature  difference  between  the 
air-preheating  tube  and  the  cell  is  relatively  constant.  At  an 
axial  distance  greater  than  ^0.34  m,  the  rate  of  radiative 
heat  transfer  decreases  slightly  from  ~700  to  ^550  Wm-2 
before  drastically  increasing  from  ^550  to  ~1450  Wm-2  in 


the  fuel  outlet  region.  The  high  radiation  rate  at  z  ~  0.5  m 
is  caused  by  the  large  temperature  difference  between  the 
cell  and  the  inlet  air  and  makes  an  important  contribution 
to  air  preheating.  Fig.  9  compares  the  cell  temperature  pro¬ 
files  in  the  cases  where  the  radiative  heat  transfer  is  included 
and  where  it  is  not  and  indicates  that  its  inclusion  leads  to 
a  decrease  in  the  cell  temperature  of  ~60  °C  near  the  cen¬ 
tre  of  the  cell.  This  temperature  difference  is  relatively  small 
because  heat  convection  in  the  air  channel  is  the  predomi¬ 
nant  mode  of  heat  transfer.  The  contribution  of  heat  radiation 
should  increase  as  the  air  velocity  is  reduced.  Although  the 
temperature  difference  seems  small  when  compared  to  the 
operating  temperature,  it  is  still  critical  for  the  heat  manage¬ 
ment  since  the  extra  amount  of  heat  generated  inside  a  cell 
can  cause  material  degradation. 

3.2.3.  Distributions  of  reversible  cell  potential, 
overpotential  and  current  density 

Fig.  10  illustrates  the  distributions  of  the  reversible  cell 
potential  (£°),  overpotentials  and  current  density  along  the 
cell  length.  E°  decreases  along  the  cell  length  because  of 
the  depletion  of  fuel  and  oxidant  downstream.  The  reversible 


(b)  z  (m) 

Fig.  11.  Surface  plots  of  distributions  of  H2  (a)and02  (b)  concentrations  for  the  base  case  model  (/g  ~  300mAcm~2,  Vceii  =  0.7  V,  Ua ir  ~  0.167  and  t/fuei  ~  0.85). 
Concentrations  are  expressed  in  terms  of  molar  fractions. 
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cell  potential  at  the  fuel  outlet  is  about  73%  of  that  at  the 
inlet  end.  This  implies  that  the  reduction  of  the  reversible 
cell  potential  is  one  of  the  major  reasons  for  decreasing  cell 
performance.  On  the  other  hand,  the  overpotentials  at  the 
cathode  and  the  anode  decrease  as  one  moves  away  from  the 
fuel  inlet.  In  general,  the  overpotential  increases  as  the  current 
density  rises.  Nonetheless,  the  overpotentials  still  decrease 
even  when  the  current  density  increases  in  the  region  from 
z  =  0  to  ~0.2m.  The  reason  is  that  the  effect  of  increasing 
current  density  is  more  than  compensated  by  the  concomitant 
increase  in  temperature  (see  Fig.  7),  which  leads  to  a  decrease 
of  overpotential.  Moreover,  the  overpotential  at  the  cathode 
side  is  greater  than  that  at  the  anode  side  over  the  entire  length 
of  cell  because  of  inherently  slower  kinetics  of  the  cathodic 
reaction  and  large  cathode  thickness  [8]. 

As  shown  in  Fig.  10(c),  the  current  density  increases 
and  reaches  a  maximum  of  ^440mAcm~2  at  an  axial 
distance  of  ~0.14m.  Thereafter,  it  continuously  decreases 
toward  the  fuel  cell  outlet.  The  minimum  current  density  is 
~50  mA  cm-2  at  the  tube  outlet.  The  key  parameters  affect¬ 
ing  the  electrochemical  reaction  rate  are  the  local  H2  con¬ 
centration  and  temperature.  At  the  fuel  cell  inlet,  high  H2 
concentration  results  in  a  high  reaction  rate.  The  current  den¬ 
sity  keeps  increasing  although  H2  concentration  decreases 
along  the  fuel  inlet  region  because  of  an  increase  in  the  local 
temperature  as  depicted  in  Fig.  7.  Beyond  0.14  m,  the 
current  density  decreases  while  the  local  temperature  is  still 
increasing  since  the  influence  of  H2  depletion  affects  the  cur¬ 
rent  density  more  strongly  than  the  local  temperature. 

3.2.4.  Distributions  of  species  concentrations 

The  distributions  of  H2  and  O2  concentrations  within 
the  cell  are  displayed  in  Fig.  11.  The  H2  concentration  is 
quite  uniform  in  the  radial  direction  over  the  entire  length 
of  the  cell  because  of  the  high  H2  diffusivity  (Fig.  11(a)). 
On  the  other  hand,  O2  concentration  at  the  reaction  site 
(cathode/electrolyte  interface)  is  much  lower  than  its  bulk 
concentration  especially  at  the  axial  position  where  the  high 
current  density  is  observed  (Fig.  11(b)).  This  is  likely  due 
to  the  fact  that  O2  diffuses  relatively  slowly  through  the 
porous  cathode.  As  a  result,  the  overall  cell  performance  of  a 
cathode- supported  tubular  SOFC  is  limited  by  the  O2  diffu¬ 
sion  through  the  thick  porous  cathode.  Therefore,  the  use  of 
a  cell  with  a  very  thin  cathode  such  as  the  recently  developed 
high  power  density  solid  oxide  fuel  cell  (HPD-SOFC)  [27] 
should  be  able  to  overcome  this  problem. 

4.  Conclusions 

A  mechanistic  model  of  a  cathode-supported  tubular 
SOFC  that  includes  charge,  momentum,  mass  and  heat  trans¬ 
port  has  been  developed.  In  this  way,  it  avoids  using  correla¬ 
tions  for  calculating  heat  and  mass  transfer  coefficients  which 
may  result  in  large  errors.  The  model  also  includes  the  effect 
of  heat  radiation  between  the  cell  and  the  air-preheating  tube. 


The  model  shows  good  agreement  with  experimental 
current-potential  curve  data  obtained  from  the  literature.  The 
distribution  plots  of  flow,  temperature,  current  density  and 
species  concentrations  contribute  to  a  better  understanding 
of  transport  phenomena  inside  the  cell.  The  model  results 
suggest  that  O2  diffusion  through  the  thick  porous  cathode 
limits  the  cell  performance  of  a  tubular- type  SOFC. 
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